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ABSTRACT 

Prediction  of  demand  and  setting  reorder  levels  by  ex- 
ponential smoothing  and  mean  absolute  deviation  is  studied 
for  an  inventory  system.   The  relation  between  the  stand- 
ard deviation  and  mean  absolute  deviation  of  demand  dis- 
tribution is  discussed.   Asymptotic  first  and  second 
moments  of  the  exponentially-smoothed  forecast  errors  (MAD) 
are  derived.   Results  obtained  from  simulation  of  several 
normal  systems  indicated  that  the  smoothing  technique  is 
inferior  to  the  classical  maximum  likelihood  estimation 
method. 
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I .   INTRODUCTION 

In  recent  years  the  technique  of  forecasting  time 
series  by  exponential  smoothing  has  received  a  considerable 
amount  of  attention  in  the  literature.   As  a  forecasting 
technique  in  inventory  theory,  it  is  particularly  suited  to 
the  situation  where  demands  are  time  dependent  and  corre- 
lated for  much  of  its  intuitive  appeal  rests  on  weighting 
the  past  in  such  a  way  that  the  most  distant  past  receives 
less  and  less  weight.   When  demands  are  not  so  related, 
smoothing  may  or  may  not  be  a  suitable  forecasting  tech- 
nique.  In  particular,  when  demands  are  independent  and 
from  a  common  distribution  such  as  for  the  so-called  fast 
moving  items,  more  classical  methods  like  maximum  likeli- 
hood may  be  preferable  because  of  the  optimal  properties 
that  are  often  demonstrable.   Even  in  such  circumstances, 
however,  exponential  smoothing  is  still  a  candidate  for  a 
forecasting  scheme  and  it  may  very  well  be  that,  while  not 
optimal  ,  the  difference  in  such  criteria  as  stockout  risk 
for  example,  may  be  negligible  compared  to  optimal  proce- 
dures from  a  practical  point  of  view.   In  that  case,  fore- 
casting rules  need  not  be  changed  for  items  differing  in 
demand  pattern  and  a  certain  uniformity  may  thus  be 
achieved. 

The  main  purpose  of  this  thesis  is  to  investigate  a 
model  of  an  inventory  system  used  by  the  Naval  Supply 


System  Command  (NAV  SUP)  following  Brown  [Ref.  l].  The 
model,  as  described  by  Ref.  1  and  Ref.  5,  is  briefly  as 
follows : 

Single  exponential  smoothing  is  an  operator  given  in 
Ref.  1  as 

t-1 

Sfc(x)  =   Oi     Y  ^   xt_k  +  ^  x0  (1-D 

k=0 

1_T_ 

where  x.  is  the  total  demand  of  the  i    period 

ol     is  the  smoothing  constant  with   0  <  a  <  1 
|S  =  1  -  a 

Intuitively  speaking  the  smoothing  operator  is  a 
method  of  forecasting  by  weighting  past  demands  in  such  a 
way  that  the  most  recent  demand  is  weighted  by  a  and  the 
remaining  demands  are  given  less  and  less  weight  as  time 
progresses. 

The  standard  deviation,  <7,  of  the  demand  distribution, 
which  is  needed  to  set  the  reorder  levels,  is  estimated  by 
the  following  three  statistics: 


e   =  X   -  S   , (x) ,   t   period  forecast  error, 


t-1 

K= a  Z  *k  K-k1  +  j8t  |eQi  (1-2) 

k=0 
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estimate  of  the  mean  absolute  deviation  of  forecast  errors, 
and 


CT  =  V  2  4  —  A  ' 

estimate  of  the  standard  deviation  of  the  demand  distribu- 
tion. 

This  estimate  of  cr,  as  explained  in  Ref.  1,  is  based 
on  the  mean  absolute  deviation  (MAD),  A,  of  the  forecast 
errors.   Reorder  levels,  based  on  the  predetermined  risk 
level  \ji    (defined  to  be  the  probability  of  running  out  of 
stock),  are  found  as  follows: 

0  =  P(X  >  \i  +   kcr) 

Here,  \x   is  the  mean  of  X,  and  k  is  a  constant  which 
can  be  determined  from  the  probability  tables  for  a  given 
risk  level,  lb.      The  reorder  level  is  then  given  by 

Reorder  Level  =  [i  +   kcr  (1-3) 

where  k  is  the  given  above. 

In  the  analysis  this  model  was  slightly  modified  be- 
cause of  the  following  considerations: 

1.   X  ,  the  demand  asked  from  the  inventory  system 

before  it  is  placed  in  operation,  is  zero.   Although  one 

may  predict  a  number  for  a  fictitious  inventory,  there  will 

be  no  observation  for  X  .   Hence  e   is  zero. 

o  o 
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2.  (x     and  $     can  be  assumed  to  be  approximately  equal 
to  zero  for  values  of  N  >  30,  if  0.01  <:  a  ^  0.99.   In  view 
of  this,  the  modified  form  of  (1-1)  and  (1-2)  used  in  the 
analysis  was 

t-1 
St(x)  =  a  £  )9k  xt_k  (1-4) 

k=0 


A   =  a  Y^   I  et-k'  (1_5) 

k=0 

In  this  model,  not  only  is  demand  forecasted  by  expo- 
nential smoothing  but  demand  variability  is  estimated  by 
means  of  exponentially  smoothed  absolute  deviations.   Such 
estimates  are  in  turn  used  to  set  safety  levels  for  inven- 
tory order  rules. 

The  model  described  above  has  been  subjected  to  many 
criticisms.   Zehna  [Ref.  1~]    has  indicated  that  the  estimate 
of  standard  deviation  might  be  very  poor  depending  upon  the 
underlying  demand  distribution.   Based  on  the  results  ob- 
tained by  Zehna  [Ref.  8],  the  analysis  in  this  paper  takes 
on  its  starting  point  on  investigation  of  the  ratio  A/<7, 
a  key  element  in  the  estimation  of  (J.   It  is  quite  evident 
that  any  weakness  in  the  estimation  procedure  for  <T  affects 
the  reorder  levels.   An  underestimate  imposes  higher  risk 
and  an  overestimate  causes  unnecessary  expenditures.   In 
particular,  the  resulting  consequences  could  be  catastroph- 
ic in  a  military  inventory  system. 
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It  is  quite  evident  that  the  probability  distributions 
of  certain  key  statistics  would  be  desirable  in  any  exten- 
tion  of  the  existing  analysis.   Although  some  of  the  dis- 
tributions required  have  been  derived,  the  distribution  of 
A  ,  a  most  important  statistic,  could  not  be  derived  ex- 
plicitly.  However,  the  first  and  second  asymptotic  moments 
of  A  were  obtained,  and  it  is  shown  that  the  asymptotic 
variance  of  A  is  finite  for  the  case  where  the  demand  has 
a  normal  distribution.   Therefore,  the  next  step  of  the 
analysis  used  simulation  techniques  to  gain  an  insight  into 
the  behavior  of  the  distribution  of  A  ,  and  to  observe  the 
effects  of  a  finite  variance.   Reorder  levels  and  their 
dispersions  along  with  the  risk  levels  attained  are  analyz- 
ed.  The  resulting  errors  due  to  the  variability  of  safety 
levels  is  examined  in  detail  and  the  consequences  are 
evaluated  from  several  points  of  view  when  demands  are 
normal.   Moreover  the  results  are  compared  to  maximum  like- 
lihood methods  establishing  that  a  considerable  difference 
exists  in  the  two  methods.   For  items  satisfying  the  demand 
assumptions  considered  here,  then,  it  is  established  quan- 
titatively that  exponential  smoothing  is  definitely  inferior 
to  more  classical  methods. 
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II.   RATIO  OF  A/(7 

The  mean  absolute  deviation  of  a  probability  distri- 
bution is  given  by 

A=  e(|  x  -  ml )  . 

For  a  normal  distribution  the  ratio  A/a  is  the  con- 
stant sj2/ft   =   0.8.   Brown  [Ref.  l]  in  his  work  on  smoothing 
and  prediction  states, 

"Hence,  the  mean  absolute  deviation  is  proportional  to 
the  standard  deviation,  and  the  ratio  may  depend  on  the 
form  of  the  distribution  p(t),  but  dependence  is  slight." 

Now,  as  indicated  by  Zehna  in  [Ref.  7],  this  is  not 
always  the  case.   For  example,  consider  a  Poisson  distri- 
bution with  parameter  X.   A/a  is  then  a  function  of  X  and 
it  may  be  much  smaller  than  0.8  for  small  values  of  X 
[Ref. 8].   However,  there  are  some  distributions  where  the 
approximation  0.8  for  the  ratio  A/a  is  a  good  one  under 
certain  conditions.   The  quote  above  is  a  rather  general 
statement  which  may  be  justified  as  follows: 

Let  X  be  a  continuous  random  variable  with  probability 
density  function  (p.d.f.)  f (x)  and  let  u  and  JL  be  the  upper 
and  lower  limits,  respectively,  of  the  random  variable  X  on 
the  real  line.   Then 


u 
A  =   E(|X    jLt|  )       |  x  -  U  |   f  (x)dx 
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or 


,u 


,u 


A=   -    f   xf(x)dx   +    U        f(x)dx   +        xf(x)dx    -    \1    f    f(x)dx. 


Hence, 


A=   2 


u 


F(x)dx 


£ 


Here,    F(x)    is    the   cumulative   distribution   of  X, 


Then 


A  =   2 

cr  "     ct 


.M 


F(x)dx 


Let      x   =   cry   +   [i.      Then, 


A 
a 


•=  =   2 


F((Ty   +  Jl)dy 
cr 


Now  suppose  that  F(x)  is  any  distribution  function  such 
that 


EM 


*^> 


1   ,   asJ-« 


where  <£(z)  is  the  cumulative  standard  normal  distribution 
function.   Then  for  large  values  of  fi/cr 

F(ay  +  (Li)  =  *(y) 
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Hence 


^S    2 
cr 

it- 

o 

j 

Y 

1 

er 

-  i  z2 

2 

e        dz  dy 


Integrating  by  parts,  the  ratio  then  becomes  A/cr  =  s/2/lTm 
There  are  some  distributions  for  which  the  above  approxima- 
tion holds;   for  example  erlang  and  chi-square.   But  for 
other  distributions  such  a  generalization  may  not  apply. 
For  instance,  let  X  have  the  following  mass  function 

0.1   if  x  =  0 
p(x)   =  ^   0.8   if   x  =  a 
,0.1   if  x  =  2a 

Then 


A  =   0.2a 
cr2  =   1.2a 


and 

A/cr   =  .182. 

The  error  made  by  assuming  A/cr  =  0.8  is  0.618,  or  77%  , 
which  is  not  negligible. 
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III.   PROBABILITY  DISTRIBUTION  OF  A 

AND  OTHER  ASSOCIATED  STATISTICS 

Of  central  interest  in  the  use  of  MAD  is  the  ability 
to  estimate  the  demand  variability.   The  variance  of  the 

forecast  error  is  always  a  function  of  the  variance  of  de- 

2 
mand,  a  .   Therefore  the  MAD  of  the  forecast  errors  is  also 

a  function  of  cr.   Since  MAD  is  a  linear  function  of  the 

standard  deviation  [Ref.  l]  and  if  MAD  of  the  forecast 

errors  is  estimated,  then  (J  can  be  estimated  as  described 

in  [Ref.  l].   But  in  this  estimation  there  are  two  points 

to  be  considered. 

1.  Invariance  Principle: 

0"  and  A  are  related  to  each  other  by  a  constant, 
i.e. , 

cr  =  cA  . 
But  this  relation  is  in  the  population,  and  it  does  not 
mean  that  the  estimator  of  cr  and  A  will  enjoy  the  same 
relationship  necessarily.   If  these  estimates  mentioned 
above  were  MLE,  then  the  invariance  principle  does  hold. 
However,  the  definition  of  A  indicates  that  this  is  not 
the  case. 

2.  The  Constant  of  Proportionality  of  A/cr : 

If  the  underlying  demand  distribution  is  not 
normal,  then  the  assumption  that  the  proportionality  con- 
stant is  approximately  0.8  needs  to  be  validated.   For  this 
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case  the  current  model  can  be  applied  by  a  slight  modifica- 
tion  on  cr,  that  is,  replacing  vff/2  by  the  reciprocal  of 

the  appropriate  proportionality  constant. 

o 

In  this  section  we  assume  that  demand  is  Normal  (a., a   ), 

where  a  is  a  constant.   Derivation  of  the  probability  dis- 
tribution  of  A  requires  probability  distribution  for  the 
following  statistics: 
A.    EXPONENTIAL  SMOOTHING  OF  PAST  OBSERVATIONS 

(Estimate  of  the  Expected  Value  of  Demand) 


t(x)  =  a  I  p  x 


t-1 

&■   v 

t-k 


k=o 


Since  all  X's  are  mutually  independent,  the 

distribution  of  S.  (x)  for  large  values  of  t  is  known  to  be 

2       Oi 

a  normal  distribution  with  mean  a  and  variance  or  (1  +  -.  .  fl) , 

±+p 

B.  FORECAST    ERRORS 

t-i-2 


k=i 


e.     .    =  X.     .    -   a     )         $K  X,    .    .     _         i    <   t 
t-i  t-i  Li  t-i-k-1 


e  •■  _  has  a  normal  distribution  with  the  follow- 


ing mean  and  variance: 

\l  =   0  for  large  values  of  t; 


V(efc_.)  =a2  (1  +  A) 


for  large  values  of  t. 
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C.    JOINT  DISTRIBUTION  OF  FORECAST  ERRORS 

Let  U  and  V  be  t-dimensional  vectors  defined 


as 


V   = 


t-1 


U     SB 


u. 


u 


t-1 


u. 


u. 


where  U, ,  ...,  U  are  constants.   Let  Y  be  defined  as 

T 
Y  =  U  V 


or 


Y  =   >  U.e. 
oil 

i=l 


Substituting  e.  in  terms  of  the  x's, 


t  i=2 

i=l  k=0 


or 


t   i-2 


I  "A  -  a  I    I  ui^  *i-,<-i 


i=l 


i=l  k=0 
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Hence, 

t-1  i-1 


ak 


Y  =  I  xt-i  C"t-i  -  a  I  r  Vi 

i=0  k=0 


Since  the  coefficients  of  the  X's  are  constants  and  the  X's 
are  mutually  independent  normal  variables,  it  follows  that 
Y  has  a  normal  distribution.   Hence  by  the  definition  of  a 
multivariate  normal  distribution  [Ref.  2],  the  forecast 
errors  have  a  multivariate  normal  distribution.   The  corre- 
lation coefficient  between  e,  .  and  e,   .,  where  i  <  j,  can 

t-i       t-j  J 

be  computed  as  follows. 
Let 


Y  =   eL    +  eL  where   i  <  j  <  t. 

t-i     t-j 


V(Y)  =  V(e.  .)  +  (V(e.   .)  +  2Cov(e.  .,  e.   .), 
v         t-i         t-j'  t-i    t-j 


or 


Cov 


<et-i'et-j>  ■  \  LV<Y>  -  v<et-i>-v<et-j>J        t3-1' 

Rewriting    in    terms    of    the   X's, 


Y   =   X.      .     -    OX.      .     ,     -    a£X  -    -    ap    x    XX 

t-i  t-i-1  t-i-2  t-j 


a^~xx      .    .  - +  x     .-ox      ._  etpx     .     . 

r  t-j-l  t-j  t-j-1  t-j-2 
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Rearranging  terms  on  the  right-hand  side, 
j-i-2 

y  =  x.  .  -  a  )   px.    .    .   ,  +  (i-oi^~:L~l)x^ 

t-i       £->  t-i-l-k  t-" 

k=0 


t-j-1 

(1,^-)  «  2,   ^t-j-l-k 
k=0 


Since  the  terms  in  the  above  expression  are  mutually  inde- 
pendent, the  variance  of  Y  for  large  values  of  t  can  be 
found  as  follows: 

v(Y)  =  a2  (i  +  a2  iz4(j"i_1)+  il-afi-^f-l   ^(1+^j:i)2 

1-jT  1  -  j8 


Subs 


tituting  V(Y),  V(e..)i  and  V(e   .)  into  (3-1) 


oveCe   ..efc   )  =  %-  (l+  A  (1-/92  ^"^^  +  (l+^-i,  2  ,  + 


■j'    2  V"   i+)8 


(l-a/9j  i  X)2  -  ^ 


|8 


Hence,  the  asymptotic  correlation  coefficient  between  e.  ■ 

and  eL    becomes 
t-D 


(1-a^  i_1)2V  1 


(3-2) 


21 


Now  consider  the  following  function, 

*■  -  ^  0 +  Oh  u-*2y~2+  <i^y)2) +  d-^-1)2)- 1 

Here   y   is    a   continuous   variable.      Then 
dfi'       «  f-  ii!l  -    „    fl    .    9fly    lorT    fl    .    9.2y 


log   jS   +   2^    log    j9  +   2/9^    log 
0 


I±if_    20<£  .        .    ,^ 


i  - 
Rearranging    the    terms, 


g-  ?  Xog^-  ^  +  f  ♦  flf  -  il#A  *  %^£  ) 


2/9z         *  *  ^p  2/9' 


or 


dp'         /9y   log/9    ^        (1+0) «  "a  Qti9y   log  ft 

dy  2  V  $         J  2/9 


Now  log  ft  is  negative;   therefore,  the  first  derivative  of 
p1,  with  respect  to  y  is  positive.   Hence  p1  ,  is  a  monotonic 
increasing  function  of  y.   Then  positive  integral  values  of 
y  also  lie  on  this  same  continuous  curve.   Actually,  p  co- 
incides with  p',  whenever  y  is  a  positive  integer.   Hence  p 
is  a  monotonic  increasing  function  of  j-i.   The  supremum  of 
p  is  found  by  letting  j-i  -*00  in  (3-2)  : 

Sup  P  =  -i±S  (l  +  ySjjj  (1  +  1)   +  i)_  1 

=  0 
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The  minimum  value  of  p  is  found  by  letting  j-i  =  1 


<>min  =  ^  (1+  IT?  CC1-1)  +  d^)2)  +  d-«)2)  "  1 


=  ^(l  +  (l+/3)a  +  j92   -  1 


=  1±A  -   i 
2 


The  minimum  value  of  j8  may  be  zero;   hence,  p    =  -0.5 

J  mm 

D.    DISTRIBUTION  OF  |  eu 

t-i 

Rewriting  |  e   .  I  as 


e,  .  =  +  7e,  . 
t-i        t-i 


The  term 
2 


e  t-i 


^(-if.) 


has,  asymptotically,  the  chi-square  distribution  with  one 
degree  of  freedom  [Ref.  3]. 
The  term 


2 


*2<1+  I7F> 
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has,  asymptotically,  the  chi  distribution.   Hence  the  mean 
value  and  variance  of  efc  j  is  given  by 


E<  'et-i  l  >  ■  CT72=s  7  f  <3-3> 


V(   I  efc_.    I    )      =  2=i  a2    (1   +  ^  )  (3-4) 


An  attempt  to  derive  the  probability  distribution  of  £  was 
made  utilizing  the  change-of-variable  technique,  and  by 
defining  the  following  transformation: 

t-1 


k=0 


Yl  =  et-l 


Y    =  e 
*t-l    1 


The  inverse  transformation  becomes 


e   =  ±  Y 
el       t-1 


e   =  ±  Y 
e2      *t-2 


t-1 
efc  -  ±(A-   ^k(±Yk) 
k=l 

This  transformation  is  obviously  not  one-to-one;   there 

are  2   sets  in  the  domain  of  the  transformation  that  are 
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transformed  into  a  single  set  in  the  image.   For  all  cases 
the  absolute  value  of  the  Jacobian  of  transformation  is 


one.   The  joint  density  of  A  ,  Y, ,  . ..,  Y   1  then  becomes 

2t  t-1 

f(y-L *t-V6)   =  Z  9(±y1,±y2,...±(6-a  X^(±Yk)) 

i=l  k=l 

where  g(x1#  x~,  .  ..,  x  )  is  the  p.d.f.  of  a  multivariate 
normal  distribution,  and  for  every  i  one  combination  of  the 
signs  is  valid.   The  distribution  of  A  then  naturally  re- 
quires the  (t-1) -fold  integration  of  the  joint  density  of 
A  ,  Y,  ,  . . . ,  Yh_-i  •   Although  every  parameter  in  the  joint 
density  is  known,  the  labor  required  to  arrive  at  an  ex- 
plicitly-defined probability  distribution  makes  the  task 
very  difficult.   Therefore,   the  analysis  concerning  the 
dispersion  of  (J  was  not  completed  because  of  the  intracti- 
bility. 
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IV.   MOMENTS  OF  A 


The  next  step  of  the  analysis,  then,  was  a  simulation 
of  the  inventory  system  in  order  to  gain  an  insight  concern- 
ing  the  distribution  of  A  .   it  was  seen  that  the  reorder 
levels  were  fluctuating  around  the  theoretical  reorder  lev- 
els, even  after  1000  observations.   This  was  related  to  the 

asymptotically  finite  variance  of  A  as  indicated  by  [Ref.7]. 

a 

Although  computer  simulations  implied  variance  of  A   does 

A 

not  vanish,  the  moments  of  A  were  derived  in  order  that  an 
analysis  of  the  simulation  results  might  be  made.   The  mo- 

A 

ments  of  A  were  derived  for  large  values  of  t  and  are, 
therefore,  asymptotic. 

A 

A.         MEAN   VALUE    OF    A  FOR   LARGE   VALUES    OF    t. 
t-1 

k=0 


t-1 

A  IT1 


»k 


E(A)    =  a     I  jB K  E(   |  et_k  |    ) 
k=0 


'  1     /  -?-      rr 

tt    N     2-CL 


B.  VARIANCE    OF    A 


t-1  t-2    t-1 


V(A)  =  a2    I  02kV(|  efc_kl  )+20i2  I       I  i?i+^Cov(|  et_.|  ,|  et_.|  ) 
k=0  i=0    j=i+l 


(4-1) 
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Cov(|et_.|,    |et_.|)    =   E(|et_.et_.|)    -   E  (|  et_.  I  )  E  (I  efc_  .  |  ) 

(4-2) 
Thus,    substituting    (3-3),     (3-4),    and    (4-2)    into    (4-1), 


MtA\  *-2     rr2         2  U  8<3  2 

v(A)   =  — —  <j    ■=—■ -    — -q  -  r  0     a     + 

t-2    t-i 

2q2      I        I   ^  +  j    E<'  et-iet-jl)  (4'3) 

i=0    j=i+l 

This  last  expression  for  V(A)  requires  the  computation  of 

E(  e,  .  e.    ).   The  following  derivations  concern 
t-i  t-j 

Ede^e^l). 

The  joint  distribution  of  e,  .  and  e,   . ,  as  indicated 
J  t-i       t-j 

in  Section  III,  is  bivariate  normal.   The  parameters  of 
this  distribution  are 


U.,    =  E(e.     .)    =   0  or. '  *  a  jl   + 


1  v    t-i7  1  V  1+/9 


^   =  E(eu    .)    =  0  cr0   =  a    /l  + 


2  v    t-j7  2  V  l+j8 

and  asymptotic  p  is  given  by  (3-2) .   For  notational  con- 
venience, let 


x  -  m     ± 


Y  =  e.   . 
t-3 
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Then  the  joint  p.d.f.  is  given  [Ref.  4]  by 


f(x,y) 


1 

e 

2 

-  °°    <    X 

-  oo     <    y 

<  oo 

<  oo 

21T 

<T   cr  v 
x  y 

/iV 

where 


1-P    *x        x  y   (Ty 


In  that  case 


00       00 

p 


If 
I  I  xyl  f(x,y)dy  dx 


_00     _00 


I  x|  Q-J  y  f(x,y)dy  +  J  y  f(x,y)dy 


dx 


o   o 

xy  f(x,y)dydx  - 

.00   "  _00 


xy  f(x,y)dydx  - 


J  ?     c 

xy  f(x,y)dydx  +      xy  f(x,y)dydx 


o   o 


However,  by  the  definition  of  correlation  coefficient, 


o   o 


E(XY)  =       xy  f(x,y)dydx  +  I      xy  f(x,y)dydx  + 


_00     _00 


00     _0  °°     „oo 

xy  f(x,y)dydx  +       x  y  f(x,y)dydx 
00  o    o 


a  a  p 
x  y 
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Thus, 


E(|  XY|  )+E(XY)  =  2 


o        o 

00 

xy   f(x,y)dydx  + 

_CQ             _00 

1 

n 

xy  f  (x,y)  dydx  J 


Making  a  change  of  variable  in  the  first  integral,  let 


z  =  -x 


u  =  -y 


then 


o   o 


00     00 


xyf(x,y)dydx  =      zu  f(-z,-u)dudz 


.00     _00 


O    O 


Since  f(x,y)  =  f(-x,-y),  it  follows  that 


P   „o 


xyf(x,y)dydx  = 


—  CO      —  CO 


zu  f(z,  u)dudz 


o   o 


Thus 


E(|xy|)  +  E(XY)  =4      xy  f(x,y)dydx 


o   o 


(4-4) 


The  right-hand  side  of  (4-4),  without  the  coefficient  4,  is 

x2  (y-b)2 


2  2     2 

2cr  2a  ^(l-pz) 

x  p        x 

xy  f (x, y)dydx  =   |  ^  ,. , ,.  •    ^ —  dydx 


o  o 


<W2*     o  ^2^  cryVl-p2' 


(4-5) 
where  b  is  the  conditional  mean  of  Y  given  X,  and  is  given 
by  [Ref.  4], 


b=^  PX 
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The  conditional  standard  deviation  of  Y  given  X  is  given  by 


c  =  0-  Jl-p' 

y 


Let 


y  -  b 


=  z  ,  then 


x 


2(Tx    " 
xe         ye 


(v-b) 
2c 


O  <V2*      O   C^ 


X 


dydx  = 


oo     or    oo 
'  xe   x 


o  "J7*  -Jb  J*8 


(cz+b) e     .   , 
N  r-—,   ' dz  dx 


Integrating  the  inner  term  first: 


(4-6) 


b' 


cz+b     2   , 
e      dz 

-b  ^2ff 

c 


J2n 


e   2c   +  b  (1 


•(-!» 


where  $(x)  is  the  cumulative  normal  distribution.   Sub- 
stituting this  expression  into  (4-6) , 


x 


2(7 


xe 


o  <J  >J2it 
x 


ye_ 


2  2 

_  (y-bV 
2c2 


o  c  72^ 


dydx 


x 


2      2 

oo      2cr      2c 
x 
ex  e 


'  o    0"  x  2tt 

2 
x 


2(7 


2(7 


bx  e 


x 


o    a  7277 

x 


dx  - 


•  b  *(-  *   )e 

c 

o  a  ^2tt 

X 


X 


dx 


dx  + 


(4-7) 
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Denote  the  terms  on  the  right-hand  side  as  1,  2,    and  3, 
respectively,  and  compute  each  one  separately  as  follows: 
1.   Substituting  values  of  c  and  b  into  1, 


x 


«2  2 
P  x 


ex  e 


2<7x    2(l-p2)CTx2      ,°°cr„yi-p2 


2*7(7 


dx 


X. 


x  e 


2(7  2(l-p2) 

A 


x 


2*ror 


x 


dx 


Y. 


vjl-p' 


2n 


*xu-p'> 


2  1 
Vx(1-p  )2 
2ff 


2.   Substituting  the  value  of  b  into  2, 
2 


x 


x 


2a 


a00    per  x     e 
"2 

O       <7       \f2TT 
X 


X 


cr   p 
dx   =  -*=- 


2cr 


2 

x      e 


x 


dx 


cr     J2ft      o 
x 


a   cr   p 

=    y  x 

2 


3.      Substituting    the   values   of  c   and  b   into   3, 

2 


x 


2(7 


# 


-    Px 

r *r> 


O  CT    Jl-p2 

X  r 


x   e 


x  (7 


(7    727T 
X 


x    p   = 


X 


(7     p 

Y 


c00    n  2CT 

2  x 

x      e 


(7       727T        O 
X 


J    „     A/2ff' 


dt   dx 


where 


P   = 


Px 


(J^/l-p 


2' 


(4-8) 
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Let  x/cr      =  y   in    (4-8)  ;       then    (4-8)    becomes 


a  a   p 
^-*-      I    y2    e 


-Y  /2 


-r 


-°°v/2ff 


dt   dY 


where 


r  = 


PY 


^CT 


Integrating  by  parts, 


2 
_  £1 

du   =  y  e  dy 


v  =  y 


-r 


-»     J2tt 


dt 


u      =   -   e 


2 

_  2_ 
2 


,-r 


r 
2 


dv  = 


-°°  y2V 


dt-y 


J^ 


dy 


Eq.     (4-8)    becomes 


K 


2  tA 

-y  e      2     J        e_      dt 
o     \j2ft 


00  00 


-K 


Pv   e 


2-s/l-p" 


°       °  yi-p2  y2»r 


dy   + 


p00  _  y_    r_r  _  if. 


e 
o 


dt    dy 


where 


K  =  -£_3£ 


a  a   p 

x   y 

2TT7 


32 


The  first  term  in  the  above  equation  is  zero;   hence, 
(4-8)  becomes 


<t  a  p2(l-p2)h 
2ff 


+  K 


H  2  I 


V2F 


dt  dy 


(4-9) 


Now  considering  the  second  term  in  (4-9) , 


2 
2 


-r  -■ 


00 


-°°  y^ 


dt  dy 


£ 


0.5  + 


-r 


J2it 


dt  )  dy 


(4-10) 

This  result,  due  to  the  fact  that  y  changes  over  positive 
numbers  and  p  is  between  -.5  and  zero,  makes  the  argument 
of  $  (•)  positive.   Hence  (4-10)  becomes 


1 
2 


2        Jl 

-  *-  p-r  -  V 

e        e      dt  dy 
o      '  o 


=  M 


Let 


t  =  V  sff 

y   =  J2^ 


then 


,   r-=p  f°°     -U    pS        __2 

M  -  I   £  +     e 
2  V  2 


o  v2u   o 


f  ->  _  *■      , 

I   e    v2  dv  du 


where 


s  = 


33 


r        -u    c  v  __^n 


1     /   w    ,       i  e  \  ,    ,  %  n   v 

M  =  ^  ,    t  +  - — ;         — ir  ■  /  (-D      — T  dv  d^ 

V2ff      o    u  o 


n   v 
n=0 


Since    the   series   is   uniformly  convergent,    it  follows    that 


— -,  fp      -u     £  >        2n+l  n   +  i 

1     / ir  1  e"     \,    nn  (-p)  _u _    , 

M  =   2   V   2    +  7^  J      T"     ^"D  p=!  2n+1  „:(2n+l)    du 

V27r      °  u        n=0  (7l-p2) 


=  |jf  +  -^    ^-u"  *=£- 


n      -u    . 
2n+l  u      e        du 


^27r      n=0  (Jl^p2) 


2n+l     2n+l  ni 


or, 


p       •  2n+l        • 


M  =  i7f  +  -^    I<-i>n   \  a-p2 


277      n=0  *       2n+l 


It  was  previously   shown    that    p  is   always  between   -.5    and  0 
Therefore, 


-    1     <    --£=B    <    1 

yi-p2 


Hence  the  series  in  the  above  expression  is  the  expansion 


f  Arctan  QjJ==~J   [Ref.  6]  . 


Vi-p2 
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Thus    Eq.     (4-10)    becomes 


1  J?  +  7F    Arctan  C  -7^=5) 

s/l-P 


sllir 


(4-11) 


Substituting    (4-11)    into    (4-9)    yields 


2  2    * 

a  a   p   (1-p  )  cr  ct   p 


x  V 


+ 


2L^i ri  /  s  +  JL 


2ff 


v/2ff 


2   V    2 


Arctan    (         P_^  )  ) 


v/2ir  '   J  1-p2 


Combining    the    results    for    1,    2,    and   3    into    (4-7),    and   sub- 
stituting  into    (4-5)    yields 


2   ^ 

J    xy   f(x,y)dy   dx   =  -%- 
o      o 


2w 


2  2    ^s 

cr   cr    p    (1-p  )  2        cr  ct    p      ,    .  -—7         . 

x  y                                xy        f    1  /  £    ,       1 

2  V    2 


2w 


^2* 


-  Arctan  (    ,        ^,)  J 


2w 


7l-p' 


Thus, 


CT     CT 


CT    CT     P  CT    CT 


E(|XY|)    =  4^-^1-pV  +^ 1/  Arctan    (-j|p )  )  - 


CT    CT     p 

x  y 


2cr   ct 
x  v 


2,h 


-P 


(l-pZ)  2   -    p     Arctan       (       M      )  )  (4-12) 

7l-P2 
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But  Q      and  (T   are  both  equal  to  (J  J   1  +  777  .   Hence, 
x      y  V      1+jS 

2ar  (1  +  .  f. )    ,     9  t  _0 

E(|XY|)  =  ±±fi_f  (1_p^)  2  _  p  Arc tan  (   H   ))   (4-13) 

ff  7l-p^ 


Therefore, 


t-2  t-1 


(4-14) 


(1+/J2)   »(1+*)2       £,  £.+1 


where   E(|e.   ,e,_, |)  is  given  by  (4-13). 

It  is  clear  that  to  arrive  at  an  expression  for  V(A) 
in  the  form  of  (4-14)  many  tedious  and  difficult  calcula- 
tions must  be  made.   In  fact,  in  order  to  arrive  at  an 

a 

explicit  value  for  V(A) ,  one  must  deal  with  an  arctangent 

2 
term.   However,  because  (7   is  common  to  each  term,  the 

ratio  V(A)/or   can  be  defined  in  terms  of  the  known  con- 
stants, and  an  approximate  numerical  answer  can  be  found 

a    2 
for  the  ratio  V(A)/cr  .   The  computer  program  in  Appendix 

B  computes  the  first  two  terms  in  (4-13)  exactly,  and  the 
third  term  approximately.   Briefly,  the  computation  of  the 
third  term  was  accomplished  as  follows: 

Eq.  (4-13)  is  a  bounded  number  for  a  given  a  and  0. 
It  cannot  be  made  greater  than  some  finite  number  d,  what- 
ever be  the  value  of  p    .   (Note  that  p  is  the  only  variable 
in  (4-13)  that  depends  on  i  and  j,  and  it  is  bounded.) 
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Therefore,  in  the  double  summation  of  (4-14) ,  for  some 
finite  value  of  i  and  j ,  higher  order  terms  will  get  smaller 
than  a  specified  number  m  which,  in  this  case,   was  chosen 
as  10    .   The  program  was  written  so  that  the  higher  terms 
smaller  than  10    were  not  included  in  the  result,  and 
were  subsequently  dropped  from  the  summation.   The  upper 
bound  of  the  error  introduced  by  omitting  the  higher  terms 
from  the  summation  can  be  found  as  follows: 
E(|XY|)  is  bounded  from  above  by 


4  it 

d  = 


ir(l+/9)   2   1+0 

Since  the  summation  is  infinite, 


00     00 


Upper  Bound  of  Error   =  2a2  -£=   10  10  £  £  &1+^ 

i=0  j=i+l 


4tt2    10" 

-10 

P 

(1-/92) 

2 

4/9 

10 

-10 

(1+0)  2 

First,  this  error  is  additive  and  hence  is  added  to 
the  result.   Second,  whatever  be  the  value  of  0  ,  it  is 
always  much  smaller  than  one.   Hence,  the  results  contained 
in  Table  1,  obtained  by  the  computer  program  for  various 
values  of  a  were  almost  exact. 
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SMOOTHING 
CONSTANT 

VARIANCE  OF  (A)  /(J 

STANDARD 
DEVIATION  OF 

(a) /or 

0.10 

0.0204 

0.1745 

0.15 

0.0325 

0.2173 

0.20 

0.0461 

0.2553 

0.25 

0.0614 

0.2905 

0.30 

0.0785 

0.3239 

0.35 

0.0979 

0.3561 

0.40 

0.1196 

0.3876 

0.45 

0.1440 

0.4187 

0.50 

0.1716 

0.4495 

0.55 

0.2026 

0.4802 

0.80 

0.4267 

0.6341 

Table  1.   Variance  of 


A 


and  Standard  Deviation  of  — 
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V.   COMPARISON  OF  EXPONENTIAL  SMOOTHING 
AND  MLE  BY  COMPUTER  SIMULATION 

For  the  reason  explained  previously,  simulation  was 
utilized  to  explore  the  parts  of  the  model  where  theoreti- 
cal work  was  very  difficult,  and  was  also  used  to  compare 
the  exponential  smoothing  with  MLE.   The  reason  that  MLE 
was  chosen  as  the  alternative  procedure  was  because  it  has 
the  property  of  an  asymptotically  vanishing  variance.   As 
indicated  in  [Ref.  8],  the  maximum  likelihood  estimator, 

A 

a    ,  of  (J  is  given  by  the  formula 


n  .        -s2 
(x. -x) 


a     = 


y  ^i' 


n 
i=l 


where  x  is  the  sample  mean.   The  statistic  a   has  a  chi- 
square  distribution  [Ref.  3];   in  fact, 


, s-i     — -z — i  r(-^-) 


and 


n§)    N  2 


Vi 


Co)   =  fl=l  a2     (l-^rC   -IV)  ) 
v      n      v     n-1  v  T^n-lv  J 

From  these  formulas, 


re-?) 


Lim    E(a)  =  a 

n  ->0° 


Lim  _mV(cr)  =  0 

p  — >  CO     *    * 
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The  computer  program  used  for  the  simulation  is  included  in 
Appendix  A.   The  simulation  procedure  and  ensuing  results 
are  as  follows: 

One  thousand  normally-distributed  random  numbers,  re- 
presenting the  demand,  were  generated  by  a  pseudorandom 
number  generator.   Various  values  of  the  parameters,  \l,    <7, 
and  i|)  were  used.   For  every  triplet  of  the  parameters,  re- 
order levels,  risk  levels  attained,  estimates  of  <7  and 
sample  mean,  and  sample  standard  deviation  of  0"  and  of  the 
estimate  of  0"  by  MAD  were  computed.   The  value  of  Oi   was  0.2 
for  all  cases.   Table  2  presents  the  results  achieved  by 
the  simulation.   As  Table  2  indicates,  the  theoretical 
values  of  the  standard  deviations  of  both  estimates  of  (7 
were  in  complete  agreement  with  the  sample  standard  devia- 
tions of  these  same  estimates.   The  variance  of  cr  is  higher 
than  that  of  cr  in  every  case.   The  last  two  columns  give 
the  mean  square  error  (MSE)  of  (7  and  (7.   The  MSE  is  given 

A 

by  the  variance  plus  the  square  of  the  bias.   Although  <7 
is  biased  for  1000  observations,  its  bias  can  be  neglected. 
Then  the  MSE  of  cf  and  (7  are  determined  by  their  variances. 
The  comparison  of  the  MSE  columns  are  quite  striking;   the 
MSE  for  cr  reaches  values  as  high  as  100  times  that  for  (7. 
By  noting  the  similarity  between  the  theoretical  values  and 
the  simulation  results,  the  effect  on  reorder  levels  was 
analyzed  as  follows: 

Computation  of  reorder  levels  was  made  by  utilizing 
Eq.  (3-3) .   First,  the  theoretical  values  of  the  reorder 
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levels  were  computed  for  every  fixed  parameter  triplet  \i, 
<J,    i/)  .   Then  the  reorder  levels,  found  by  smoothing  and 
MLE,  were  compared  with  the  theoretical  values.   This  com- 
parison is  summarized  in  Table  3.   The  value  Q,  in  Table  4 
is  the  percentage  (rounded  off  to  the  nearest  integer)  of 
1000  reorder  levels  found  to  be  within  one  unit  of  the 
theoretical  reorder  level.   The  first  column  under  each  of 
the  headings  Q,  ,  P,-n/  and  Pq[-  lists  the  results  for  smooth- 
ed estimates;   the  second  column  under  each  heading  lists 
MLE.   The  quantity  P_n  represents  the  number  of  units  about 
the  theoretical  reorder  level  within  which  50  per  cent  of 
the  computed  reorder  levels  were  found;   P  _  is  a  similar 
computation  for  95  per  cent  of  the  reorder  levels.   Table 
3  clearly  illustrates  the  superiority  of  MLE  compared  to 
exponential  smoothing.   For  the  case  where  smoothing  was 
used  to  estimate  cr,  reorder  levels  fluctuated  around  the 
theoretical  reorder  level,  causing  either  an  excessive 
purchase,  or  imposing  more  risk  than  desired.   The  magni- 
tudes of  these  fluctuations  are  clearly  evident  from  the 
table  . 

Table  4  indicates  the  risk  levels  attained  by  both 
methods  at  the  1000th  observation.   For  every  predetermined 
parameter  pair  \x   and  cr,  five  different  risk  levels  were 
selected.   For  every  triplet,  the  theoretical  reorder 
levels  were  compared  with  the  computed  reorder  levels  of 
both  procedures.   The  actual  risk  levels  attained  were  then 
computed. 
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PARAMETER 
TRIPLETS 

Ql 

P50 

P95 

40,4, .01 

27 

95 

2 

0 

6 

0 

40,4, .11 

41 

94 

2 

0 

4 

0 

40,4, .50 

49 

92 

1 

0 

3 

1 

50,5, .01 

22 

96 

3 

0 

8 

0 

50,5, .11 

32 

98 

2 

0 

5 

0 

50,5, .50 

46 

99 

1 

0 

3 

0 

100.10, .01 

12 

90 

5 

0 

14 

1 

100,10, .11 

18 

94 

4 

0 

11 

0 

100,10, .50 

26 

94 

2 

0 

6 

0 

200,20, .01 

5 

71 

10 

1 

29 

6 

200,20, .11 

9 

82 

6 

1 

18 

3 

200,20, .50 

10 

84 

5 

1 

14 

2 

300,30, .01 

4 

81 

12 

1 

37 

6 

300,30, .11 

7 

74 

8 

1 

24 

3 

300,30, .50 

8 

51 

7 

1 

18 

3 

400,40, .01 

2 

15 

20 

3 

53 

10 

400,40, .11 

4 

32 

13 

2 

37 

6 

400,40, .50 

5 

75 

10 

1 

27 

3 

500,50, .01 

3 

1 

24 

9 

58 

11 

500,50, .11 

3 

4 

15 

4 

38 

6 

500,50, .50 

4 

73 

11 

1 

29 

4 

600,60, .01 

1 

12 

27 

3 

77 

17 

600,60, .11 

3 

10 

18 

3 

51 

5 

600,60, .50 

3 

1 

14 

4 

42 

12 

700,70, .01 

1 

33 

32 

3 

93 

16 

700,70, .11 

3 

35 

20 

2 

63 

10 

700,70, .50 

4 

16 

15 

2 

45 

15 

800,80, .01 

1 

0 

37 

12 

98 

27 

800,80, .11 

2 

0 

25 

8 

64 

18 

800,80, .50 

3 

32 

18 

2 

54 

9 

Table  3.   Percentiles  for  Reorder 
Levels. 
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PARAMETER 
TRIPLETS 

Ql 

P50 

P95 

900,90, .01 

1 

30 

42 

2 

125 

12 

900,90, .11 

2 

37 

29 

2 

83 

12 

1000,100, .01 

1 

25 

45 

3 

126 

12 

1000,100, .11 

1 

8 

29 

4 

87 

7 

1000,100, .50 
1 

1 

6 

25 

2 

68 

13 

Table  3 . (Continued)  - 


Percentiles  for  Reorder 
Levels. 
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As  Table  4  indicates,  the  risk  levels  attained  by  MLE 
are  found  to  be  very  close  to  the  desired  risk  levels. 
However,   the  risk  levels  of  the  smoothing  method  were 
found  to  fluctuate  over  a  wide  range. 

These  results,  gained  by  the  simulation,  indicated 
that  the  MLE  technique  is  to  be  preferred  over  the  tech- 
nique of  exponential  smoothing.   However,  due  to  the  lack 
of  theoretical  analysis  concerning  the  distribution  of 
these  results  are  not  considered  conclusive.  However,  this 
analysis  does  indicate  that  whenever  a  constant  mean  de- 
mand can  be  justified,  the  MLE  technique  is  to  be  prefer- 
red to  the  exponential  smoothing  technique. 
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APPENDIX   A 

ST  MULATION 

C  IN    THF    FALLOWING     PROGRAM    TWO    M!=THQDS    OF    ESTIMATION    OF 

C  MFAN    ANn    STANDARD    DEVIATION    0<=     A    STATIONARY    PROCESS 

f  *RE    CONT&ASTFD    BY    USING    SIMULATION    TECHN?QUFS 

C  METHOD*    ARE: 

r  l.THF    EXPONENTIALLY    SMOOTHED    D-Tfi    AS    THE    ESTIMATE    OF 

C  MEAN    AND    THE    EXPONENTIALLY    SMCPTHEO    ABSOLUTE    VALUES 

C  OF    FPROPS    AS    THE    ESTIMATE    OF    STANDARD    DEVIATION    TIMFS 

C  A     KNOWN    CONSTANT 

C  2. MAXIMijv    LIKFLIHOOD    ESTIMATES 

r  ALSO     IN    THE    PROGRAM    THE    EFFECTS     OF    ROTH    ESTIMATION 

C  PRPCEDUPFS    ON    THP    RISK    AND    «>P  ORDER    LEVFLS     ARP    <^TMIJLA- 

C  TED    BY    SAMPLE     SI  7E    1COO 

DIMENSION    "K(6)  ,T<6)  ,Xl(6)  ,XP(6),XM(6)XN(6),XK(6) 
IXlf  POO, 5)  ,  X2U0C0,5)  .G1U0Q5)  ,G2UQ05),  mO),TX<6) 
2"TY(6)  ,RISKX(6)  ,RISKYC6) 

C  RISK    LEVEL* 

PK< 1)=0.01 

«=K(  2)=C.C5 
PK( 3)=C.ll 
DK(  *)=C25 
c*( 5)=C.CC 

READ(5,131)     (T(I)fI=l,K) 
131    F0RMAT(5F4.2) 
S=10.0 
XJ=C.C 
Y=URN( 01 
DO    »21    MT=l,16 
XJ=XJ+1.C 
DO     1     J  =  l,5 
DO     1     1=1,800 
Xl<  I, J)  =0.0 
1    X2U,J)=C.C 
?L=S 

C  SMOOTHING    CONSTANT    ALPHA    I*     3.2 

BETA=O.P 

WR|TB(6i879J     S,STDFX 
97C    cqpma  T(  '1  ,TA5,  «MEAM='  ,«=6.1  ,l")Xt  «*.PCV=»  ,F6.1) 
QELTA=0.798*SQRT{2.0/( l.C+BETAJ  )**TDEX 

C  COMMUTATION    OF     THE    THEORETICAL    REORDER    LEVELS 

DO    11    1=1,5 
11    71  I  )=    7L*STDEX*TU  ) 

WRTTEl 6,150)     (7(1)  ,1=1,5) 
150    FORMAT! //T10, 'THEORETICAL    REORDER    LFVPt S • , 5F  10. 3// ) 

V=1.C 

SUM=7L 

TOT=o.C 

R 1=0.0 

*2=0.C 

B3=C.C 

B4=C.O 

C  GENERATION    OF    NORMALLY    DISTIBUTED    NUMBERS    WHICH    HAVE 

C  MEAN    S    AND    STANDARD    DEVIATION    STDEX 
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DO    21    I=l,10CO 

CAU     NORMAL    (C.0,STPFX,X) 

Y=7L+X 


METHOD     1 


OELTA=0.2*ABS(  S-  Y)  +  C.  8*DELTA 

S=C.2*Y+C.q*S 

<=DEV=l.ieP*DELT* 

V=V*1.G 

Sl!M  =  SL,M  +  Y 


METHOD    2 


X3/R=SL'M/V 

TOT=TOT+( XB&o-Y)**2 

STDEV=SOPT(TOT/( V-l.C) ) 

OK  I  )=SPEV 

G2( I  )=STDEV 

DO    20    J=l,5 

XL( J)=S+SDEV*T( J) 

XP(  .1)=XRAR-»-STDE  V*T(  J) 

XK( J)=Xl( J)-Z( J) 

TX( J)=( XL( J)-Zt) /STDEX 


RISKS    ATTAINED    BY    BOTH     METHODS 


20 


101 

147 

21 


R  T  SK  X  (  J  )  =  C  .  5-0 .  5  *E  RF  (  T X  (  J )  *0  .  7C  7  ) 
TY( J)  =  ( XD( J)-ZL) /STDEX 
RISKY(J)=C.5-0.5*ERF(TY(J) *0.7G7) 

xm(  n  =  xp(  j)-z(  j) 

M1=XK( J) 
N2=XM( J) 

IP(N1.LT.0)     M=-N1 
IFCN2.LT. 0)    N2=-N2 

XHNl  +  lffJ)-XKNl+lf  J)*1,0 

X2f N2+1,J)=X2(N2+1  ,  J)  4-1.0 

NS=I 

IF(N?.LF.899)    00    TC    21 

WR!TF(6,1C1)     CXL(M),M=1 

FORMAT(  /T13,5F7.1 f10Xt5F7.1) 

WRITF(6,147)     (RISKX(M)  ,M=1,5)  ,(PISKY(M)  ,m=i 

F0RMAT(/T2,  'ACTUAL    R  T  SK  '  ,5  F7  .  3  , 10X  ,5  F7  .  3  ) 

CONTINUE 


5)  ,(XP(V) .m=i,5) 


5) 


C  COMPUTATION    OF    THE    OBSERVED    CRFQUENCIE*    WITHIN    THE 

C  NEIGHBORHOOD    OF    ACTUAL    REORDER    LEVEL    WHEPE    1=1, 2,3, 


10 


25 

24 

35 
112 


DO 

WRI 

FOR 

Al= 

A  2= 

1  =  1 

IF( 

Al  = 

IF( 

A2= 

WRI 

cnp 
1    *A 

1  =  1 
XA  = 


26    J=l,5 

TF(6,1Q5)    PK(J) 

MAT(//T45,'PISK=«,F5.2) 

0.0 
CO 


Al.GE 
Xl(  I, 
A2.GE 
X2(  TT 

TEC  6, 

MAT(T 

NO-'  , 

4-1 

A1+A2 


.  1.0)    00    TO    24 
J)/1COO.O+A1 
.1.0)    00    TO    35 
J)/1C00.0«-A2 
112)     I,A1,A2 
10,  'PERCENTILE 
I4f2F10.4) 


BETWEEN    THEO.    REORDnR    LEVEL 
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IF(  X*  .LF.  l.<5  96)  GOTO  25 

26  CONTINUE 

C      SAMPLE  *EAN  AND  STANDARD  DEVIATION  OF  BOTH  FSTIMATES 
C      OF  POPULATION  STANDARD  DEVIATION 
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2S 


200 


2C1 


721 

722 
821 


DO 
81  = 
B2= 
3  1  = 
B2= 
DO 
B3= 
84= 
B3= 
B4= 
WRI 
POR 
IDE 
WPI 
FOP 
IMS 
S=l 
IF( 

s=x 

00 
c  =  «: 

CON 
sjo 


2^  1  =  1 
G 1  (  I  )  + 
82  +  G2< 
R1/1C0 
B2/1CC 
28  1  =  1 
B3+fGl 
R4-MB2 
SORT(R 
SORT(R 
TEC  6,2 
MA  T  (  /  / 
V  ) 

TE( £,2 
WAT(// 
2) *  ,FX 

oc.c 

XJ-5.C 
J*1C0 

TO  821 
*(  XJ-5 

TINUE 
P 


,1CC0 

Bl 

I) 

0.0 

CO 

,1000 

<I)-B1)**2 

-G2U I >**2 

3/1G00.0) 

4/1COO.0) 

00) 

TA5, 'SAMPLE     MEAN    AND    S.    DEV    OF    ESTIMATED    S 

CD    B1»B3,B2.B4 
T4*,  MMAD)  " 


) 


)    B1»B3,B2.B4 

•>,  MMAP)  »,5X,  «MFAN='  ,F7  .  2  ,  •  S  DEV=  •  ,  «=7.2,15X, 

MFAN=« ,F7.2, «SDEV=» ,F7.2) 

721,721,722 


•  01 


SUBROUTINE    NORMAL    (EX,STOX,X) 

SUM=G.C 

DO    5    1=1,26 

R=URN(  1) 

SUM=SUM+P 

SUM=SUM/26.0 

X=ST0X*17.61*(SUM-G.5)-»-EX 

RETURN 

PNO 
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APPENDIX    B 

NUMEO.TCAI     COMPUTATION    OF    VA»MNCF    OF  A 

IMPLICIT    q-*l*«    (A-H,C-Z) 
no    20    K=l,ll 

SMOOTHING    CONSTANT 

Al_PHA  =  l.-BETA 
A  ^l!M  =  O.C 

WRTTM  6,110)     *l.pHA 
110    FORMATC //T10, 'SMOOTHING    CONSTANT    = » , f^ . 3  ) 

fipcp.r  xIM.«  TTON    OF    THF     INFINITE    SUM 

no    1    M»l ,300 

I=M-1 
J=M 

T=^  ,0/(  3.1M*(1.0  +  BFTA>  ) 
P=AlDHA/( l.+BExA ) 

n=RPT<\**(  i-kj) 

C=BET4**f  J- 1  ) 
A=C/BETA 

RHO  =  <  l.-t-BETA  )/-.*{  l.  +  P*(l.  0-A**2+(  1  .+B)**2)  +  (  1.- ALPHA* 
1A  )**2-l. 

7M=P<:qR'r(  i.-PHO**2l 

7P=_PHQ/7M 

<:U^  =  T*(  ?N-RH0*r-ATAN(  7  P  )  )  *n 

I«=(  SlJW-l.OP-10)     1,1,2 
2    fi=A*BcTA 
D=0*BETA 
00   Tn    3 

1    CONTINUE 

W2TTF(A     102)     SUM 
1C2    P0PMAT(T1C,'LAST    SUM=« ,F20.12) 

X=ALCHA/(  l.+RETA 1**2*2.2  82  /3.1<-l 
7_=8.C/3.14i*BETA/( i.O+BETA)**2 

VARIANCE    Oc  A 

VAR=X-Z+2.Q*(1. O-BFTA ) **2*ASUM 

a/MPA^    O^VIA'ION    np  ^ 

?DEV=OSOPT( 3.141 *(1.+BFTA )*VAR) /2.C 

WP!TE(6,1GD     VAR,SDFV 
1C1    FQRMAT(T10,  'VARIANCE    OF     VAO • , Fl 5 .7 , 5 X , ' ST AND ART     DEVIA 
1TION    SF15.7I 
RETA=BFTA-.05 
IF(K.?0.1C)     BETA=.2 
20    CONTINUF 
cjOP 
END 
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